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ABSTRACT 

We present an attempt to improve models of the Rossiter-McLaughlin effect by relax¬ 
ing several restrictive assumptions. We consider the entire multiline stellar spectrum 
rather than just a single line, use no assumptions about the shape of the lines pro¬ 
files, and allow arbitrary size ratio for the star and its eclipser. However, we neglect 
the effect of macro-turbulence and differential rotation. We construct our model as 
a power series in the stellar rotation velocity, Vsin/, giving a closed set of analytic 
formulae for up to three terms, and assuming quadratic limb-darkening law. We con¬ 
sider three major approaches of determining the Doppler shift: cross-correlation with 
a predefined template, cross-correlation with an out-of-transit stellar spectrum, and 
parametric modelling of the spectrum. 

A numerical testcase revels that our model preserves good accuracy for the rotation 
velocity of up to the limit of 2 — 3 times the average linewidth in the spectrum. We 
also apply our approach to the Doppler data of HD 189733, for which we obtain an 
improved model of the Rossiter-McLaughlin effect with two correction terms, and 
derive a reduced value for V sin/. 

Key words: techniques: radial velocities - methods: data analysis - methods: ana¬ 
lytical - planetary systems - stars: individual: HD 189733 


1 INTRODUCTION 


Whereas the number of the discovered exoplanets grows continuously, the importance of their cross-characterization by inde¬ 
pendent observation techniques increases. There are two mostly productive planet detection methods: by radial velocity (RV) 
variations and by a photometric fadening during a transit. Consequently, the joint analysis of the combined RV-|-transit data 
gained a special value in the recent years. 

This task is not reduced to a mere combination of the RV and transit data, with their respective separate models. In such 
cases we may also observe hybrid events, like the Rossiter-McLaughlin (RM) effect, which is basically a spectrosopic view of 
a planetary transit before a rotating star. 

The most simple model of this eff ect is based on the assumption that the mea sured Doppler anomaly is equal to the 
average RV of the occulted stellar disk I Kopall 19421 : Ohta et al. 200!tI : Gimenez 200fil) . We will call this as classic RM model. 
This average velocity is given by an exact formula 
fVp 


Hnean — 


1 -/’ 


( 1 . 1 ) 


where / is a fraction of the flux blocked by the planet, and Vp is the average “subplanet” RV, computed with an account for 
the stellar limb darkening. The remaining problem here is to compute / and Vp for an assumed limb-darkening law. However, 
Hnean In this formula is not the same physical quantity as the Doppler anomaly that we seek. Instead of averaging the RV 
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2 R.V. Baluev & V.Sh. Shaidulin 


over the stellar disk, we must average the stellar spectrum first, and then determine the Doppler anomaly from this average 
spectrum. 

The average of a function is not equal to the function of an averaged argument, unless the function is linear or well- 
linearizable. Therefore, the formulae HU approximates the RM anomaly only if the stellar rotation velocity Vsini is very 
small. Ideally, it should be much smaller than the typical width of the spectral lines (scaled in the velocity units). In this 
case, stellar spectrum can be linearized with respect to the rotational Doppler shift. In the remaining cases, the formula inj 
cannot be used to predict the RM anomaly, with enough accuracy at least. _ 

There are works in which an attempt is made to construct more accurate approximation than inj. see e.g. dHirano et al.l 
20inl . 2011 : Boue et al. 201.lll . They succeeded a lot in this field, but the problem is still far from being solved due to some 


restrictive asumptions adopted by these authors. A major one is that their analytic results refer to a simplified single-line 
model of the stellar spectra. In practice, however, Doppler shift is determined from rich spectra that contain thousands of 
lines or more. Assumption of a single line cannot leav e no implications on the reliability of the model. Another important 
assumptions are that line profiles should be symmetric (IBoue et al.ll2013h and the planet is assumed small in all these works. 

Our aim here is to consider the full stellar spectrum containing multiple lines, and discuss the differences with a single-line 
model. Also, we tried to avoid decompositions in planet radius, whenever possible. This may be useful for red dwarfs transited 
by a giant planet. In this case, the planet/star radii ratio may exceed 1/10. As far as we could l earn, the largest or one of the 
largest values for this ratio currently belongs to the unique circumbinary planet KIC 9632895 ( Welsh et al. 201511 . Here this 
ratio reaches 0.26 for one of the binary components, although the absolute planet radius is only ().2Rq. In theory, a red dwarf 
star can be even smaller than a giant planet, so such a ratio can even be comparable to or exceed unit. 

The structure of the paper is summarized as follows. In Sect. [2] we give general mathematical formulation and present 
several methods and main formulae that are useful for analytic modelling the RM effect under different assumptions. In Sect. [3] 
we derive our main results of the RM effect. In Sect.|4]we describe an analytic computation of the RV momenta in an occulted 
stellar disk, which appear in our RM model. In Sect. [5] we present results of a simulation to test the accuracy and usefullness 
of the model. In Sect. [6] we apply our models of the RM effect to the public data of the MS star HD 189733, using it as a 
testcase. 


2 MAIN MATHEMATICAL METHODS AND TECHNIQUES FOR MODELLING THE 
ROSSITER-MCLAUGHLIN ANOMALY 

2.1 General formulae and definitions 

Let us adopt the logarithmic scale in the wavelength, i- = InA, and denote the spectrum of the star surface near the disk centre 
as The Doppler-shifted spectrum should then be f), where the Doppler shift is f = Vj/c-|- jB), with Vj being 

the radial velocity of an emitting point (the z axis directed along the line of sight). This is a non-relativistic approximation. 
Light coming from different points in the visible stellar disk is combined with different Doppler shifts and different local 
brightness, forming two auxiliary spectra: the cumulative star spectrum and the “subplanet” spectrum which is 

generated by a portion of the surface blocked by the transiting object. These spectra can be expressed as follows: 

= j ^{s — vx)I{\R\^s)dR, = J ,^{s—Vx)I{\R\,s)dR, R = {x,y}. (2.1) 

where t) = is a renormalized rotation velocity, I{R,s) is the limb-darkening law normalized to 7(0) = 1. This law may 
depend on the wavelength. The integration is done either over the entire star disk |i?| < 1 or over the subplanet portion of the 
disk o9p. The star radius is assumed unit here, meaning that radial velocity of each point of the surface is equal to just vx. 
The observed star spectrum during a transit is then expressed as ^t{s) = — The formulae 1231 ) assume that their 

integrad does not depend on the point in the stellar disk, except for via the rotational Doppler shift and limb darkening law 
that may vary with wavelength. Some effects may induce additional changes. For example, macro-turbulence in the stellar 
atmosphere makes lines ch aracteristics different in the disk cen tre and near the limb due to different projected geometry of 
the turbulent motions (e.g. Hirano et al. 2011 : Boue et al. 2013h . Here we do not take into account effects of this type. 

From im, the spectrum of a non-rotating star would be 


^^{s)=^{s) J I{\R\,s)dR, 


( 2 . 2 ) 


|Ji|<l 


which slightly differs from the surface spectrum ^ ( 5 ) due to the wavelength dependence of the limb-darkening law. However, 
the multiplier near is a slowly varying function, so in practice the difference between and is not important. Below, 
we will often say “non-rotating sta r spectrum” actually meaning 

Contrary to Boue et al. ( 20131 1. we do not make an assumption that ^{s) contains only a single line, and also we honor 
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Analytic models for the Rossiter-McLaughlin effect 3 

the dependence of the limb-darkening law on the spectral range. Concerning the notations, we do not introduce an explicit 
Doppler shift to the argument of at this stage, and we do not normalize onr spectra to unit. 


2.2 Modelling the procedure of determining the Doppler shift from the spectrum 


Now assume that we have a comparison, or template, spectrum and seek the best fitting Doppler shift by minimizing 

the goodness-of-fit function as follows: 


-poo 

X^{s,a) = j {^t{s)-a^T:{s-s)fdsv 


■ mm. 

s,a 


(2.3) 


in the sense of the L 2 metric. With 


From now on, let us introduce the scalar product of functions {f,g) and the norm | 
these definitions we may write the following: 

X^{s,a) = (2.4) 

The first and the third terms here do not depend on s, so to fit s means to maximize the cross-correlation function (CCF): 

—s)) \—>■ max => — s)) =Q. (2.5) 

S 

Note that without loss of generality we may assume that CCF of with is maximized at f = 0, implying that 

= 0 . ( 2 . 6 ) 

This means that the template is centred so that for an uneclipsed star the fitted RV is zero, and during the transit we 
deal with only the RV offset due to the RM effect Q 

The comparison template may be either an a priori given mask, or we may adopt it be equal to the uneclipsed star 
spectrum Up to a certain degree, these cases model the classic CCF a pproach, an d to the iod i ne cel l technique of 

Doppler measurements, respectively ( Hirano et al. 2nifll . 20 111 : Bone et al. 2013h . Note that Boue et al. ll201,3ll say that the 
line profile should be symmetric to have (lOll be equivalent to (I231l . We believe this requirement is excessive, because their 
integral I 2 in their eqs. (18,19) is always zero, even when the profile is asymmetric. This can be established by integrating it 
by parts, taking into account that boundary effects are negligible if the total integration range is large (they actually assume 
it is infinite). In fact, this integral is equal to the derivative of our ||,^t||^ over s, but this norm does not depend on any shift 
in the integration variable (again, neglecting the boundary effects). 

Of course, the practical procedures of determining s are always more complicated than in the approximations adopted 
above. For example, when dealing with a predefined template mask the CCF CtT(i) is actually not directly maximized 
but first fitted by a Gaussian '^csis — s) via s and (7, and then the fitted value of s is adopted as a Doppler shift estimate 
i Baranne et al. 19961 : Pepe et al. 2 OO 2 I '). In this case we should solve a secondary ;|j^-minimization task: 


> mm. 

s,(y,a 


XcCFi^y'^y'^) = j {CtRs) - a'^aL - s)f ds = \\Ct'i\f-2a{CMs)<:^a{s - S)) + af\'^a\\ 

Obviously, finding s is again equivalent to maximizing just a CCF, but already of a second-level one: 

-|-oo 4-®° 

CtT('f|o') = (CtT('5)^(T(7 — 7 )) = JJ—— s)duds = —s\(y)) , = J 


(2.7) 


'/jdu. (2.8) 


As we can see, this method becomes equivalent to the one with direct CCF maximization, if we replace the original template 
by convolved with the fitted Gaussian (thus imposing some broadening effect on the lines of But now it becomes 
important that depends on the parameter a, which should be fitted simultaneously with s. The best fitting values of a 
and <T can be obtained by equating the partial derivatives of (E3 to zero. Taking into account the transformation (iTSll we 
finally obtain an implicit equation for (T: 


- s\(y)) = ■ 


(.^,(.)#{'(.-f|(T)) = - 


2ct‘^ = -- 


(2.9) 


2(7v^’ X'" 1 ‘■1”^/ 4cT3yi (.^t(7)#x(7 -j|ct)) ' 

Here we used an identity , implying that d^-Yjdo = By a convention, the stroke always refers to 

derivatives with respect to 7 , not a. Additionally, instead of (ESI), we must satisfy analogous equation for .’Pi'. 

(.^*(7)#t(7|cto)) 


(.^*(7)#{(5|(To))=0, 


2ol = 


( 2 . 10 ) 


(£4f)^^ko)) 

The iodine cell techniques iButler et al. 19961 : Anglada-Escude fc Butler 2012 ') are also much more complicated then 


^ In either case, we do not take into account the Doppler shift due to the motion of the star around the star-planet system barycentre. 
We always consider only Doppler shifts relatively to the star orbital motion. 
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4 R.V. Baluev & V.Sh. Shaidulin 


the simplified fitting like (12.311 . In particular, the number of spectral parameters is much larger than two. In this work we 
adopt (1^ with ,'^1 = as an approximation to the reality. In this approximation, the resulting Doppler shift should become 
the same as if we plainly maximized the CCF with The method of cross-correlating with a reference sta r spectrum is, 
by the way, another independent Doppler technique that is used in practice sometimes (ILanotte et al.ll2014n . 


2.3 Two types of approximations leading to a “small” RM anomaly 


To move any further from (1^ . we may need to assume that s is small enough to justify the power series decomposition in s. 
This assumption becomes valid when one of the following is satisfied: 

(i) Rotation velocity V sini is small enough in comparison with typical line widths (in the spectrum of a non-rotating star). 
In this case we can decompose both spectra and ,^p into powers of t). Regardless of this restriction, the size of the 
transiting object can be arbitrary here, e.g. comparable to the star itself or even larger. Also, this approach does not need to 
make assumptions about shapes of spectral lines. 

(ii) Relative flux drop / during the transit is small enough, so that causes only a small anomaly to each line in the 
combined spectrum However, this anomaly may be shifted significantly, even by a quantity larger or much larger than 
typical line widths for a non-rotating star. In this case we can introduce various power-series decompositions in /, but we 
cannot decompose and hence So, we have to either use more or less realistic approximations of t he line shapes (e.g . 
assume they are Gaussian) or to use numeric computations where required. This is the approach adopted bv iBoue et al.l 1 201311 . 
Note that in this method it is still legal to decompose J^p into powers of Doppler shift after a proper centering, because if the 
planet is small it blocks only a small range of surface rotation velocities, well below the typical line widths. 

Note that e.g. Hirano et al. ( 201Cll l uses both these assumptions simultaneously. 

Regardless of which of the above assumptions is adopted, let us first handle the necessary decomposition of in dlSl): 

- f(^t^T) + J ^ = 0- (2.11) 

Note that by using (1^ we may derive that Also, we may perform an integration by parts in any scalar 

product of the type to move differentiations from one its operand to another, when necessary. 

The solution for s can be derived from (ItTTIi by successive approximations, and the first three terms look like: 

„ (^p^ip) 1 / {^p^!p) \" 1 / {^p^{) \ ^ f ) y 1 / .2121 

■ ^ 2 

The first-order approximation j is a small quanitity, so (12.1211 represents actually a power series in fi. 

Its error is then ^(^)- 

If the Doppler shift is determined by fitting the CCF with a Gaussian, as in (123), we should replace with ^-p and 
also need to provide an approximation for two variables f. O'. To reach this goal, we consider the system of two equations, 
Cjj(i') = 0 and the last one in (12.91) for o, and linearize them about the point f = 0 and o = Oq. Taking into account (12.101) . 
this yielded the following first-order approximation: 


((.^p#x) +2o2 (.^P#")) 






2o2 








(2.13) 


For shortness, without arguments corresponds to O = Oq here. In what follows below, we do not need more terms in the 
decomposition (I2.13t . 


2.4 Comparison with ( Boue et al. 2013h 

Bone et al. I (l201.3h assume that ^p oc f and use only the first-order approximation in /. In this case our formula (l2T2ll can be 


reduced as follows 




and formulae (l2T^ turns in the similar way into 

(j2 - (72 (,^p#x) - {■^P'^t) 






(2.14) 


(2.15) 


2 5 - 2ct2 (.^*#7 

Furthermore, iBoue et al.l (l2ni,3h consider the model with only single-lined spectra, in particular a plain Gaussian profile ir 
Also, they consider that line profiles are symmetric. This necessitates that = 0 and hence s is finally expressed 

by almost the same formulae in (12.141) and (12.1511 . while the value of CT becomes not important. 
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Analytic models for the Rossiter-McLaughlin effect 5 


However, Boue et al. ( 2013l l do not mimic the procedures of iBaranne et all ( 1996h : Pene et ^ ll2002li strictly. Instead 
of constructing the CCF with a predefined template and subseque nt fit of this CCF by a Gaussian, t hey a ssume that the 
template is a fittable Gaussian itself. We follow the sequence by Baranne et al. ( 1996li : Pene et al. ( 20021') more strictly, 
considering no fittable parameters in but instead performing a Gaussian fit of the resulting GGF. Therefore, our results 
for the GGF technique should not n ecessarily coincide with those by Boue et al. ( 2013l i in general. Nonetheless, it is possible 
to bridge them. Using eqs. (6) from iBoue et al. ( 201.ll i for the case f = 0 (out-of-transit state) we can derive that = 

= — aulA(R.\fTl). In this formulae, ag and (Jq are the best-fit parameters of the Gaussian template, as 
, and is differ ent from our definition. This additional relation allows us to reproduce entirely the 


= -aoII^tI 

2 = -C 

Boue et al. 

20131 


n formula (12) from ( Boue et al work, based on our formula arm . 

Whenever coincides with we obtain from (12.141) : 


(2.16) 


Taking into account all differences in the notation, this replicates eq. (27) from ( Boue et al. 2ni3h that expresses the RM 

anomaly for the iodine cell technique. _ 

Thus, our formulae allow to confidently reproduce the main results from ( Boue et al. 2013li . but rely on more general 
formulations (except for the effect of macro-turbulence that we neglect). 


2.5 Approximations of the star and subplanet spectra 


Let us first provide the least restrictive decomposition for 
^p{s) = j ^{s-vx)l{\R\,s)dR = j 




^{s-Sp) + ^’{s-Sp)(Sp-Vx) +-^"{s-Sp){Sp- 


vxf 


I{\R\,s)dR = 


1 


^{s-Sp)Mo+^'{s- ip) (ipMo - vMi ) + Sp)(v^M2 - IvspMi + SpMg) - 


Mi^(s) = j Al(\R\,s)dR. 


(2.17) 


Note that M]^ defined above are unrelated to Mj. defined by Boue et al. ( 2013li . Gurrently the decomposition point ip is rather 
arbitrary, and we still can choose it as we like. We may notice that whenever ip = vM\/Mq, the linear term in this series 
vanishes, leaving only the quadratic and higher te rms. Therefore, th is is the natural reference point for the decomposition. It 
coincides with the “subplanet velocity” defined in ( Boue et al. 2013h . Thus, we can write down 


Mo(i) 


= .^(i-ip(i)) + -.^"(i-ip(i))fTp^(i)--.^'"(i-ip(i))rp(.) + ^(u\0. 


,(i) = V 


Mi(i) 

Mo(i) 


_ ,,2-^2(‘S( 2 


oRs) = 




rp(s) = 


(2.18) 


This decomposition of ^p remains equally valid for the both limiting cases introduced above, slow rotation or small planet. 


In both these cases, CTp and 7p appear small. The argument Ux —i-p in (12.171) is either of the order ff{v) or 


(r), where r is the 


planet/star radii ratio. This implies that CTp is either ff{vA or and 7p is either or and the remaining terms 

of (12.181) are by an order higher. The quantity i-p is either ^(v) or so it is small only in the case of slow rotation. In fact, 

the first two terms in (I2dll) only reflect the effect of Doppler shift by ip and the rotational line broadening effect, characterized 
by CTp. The third term characterizes the asymmetry effect of the rotational broadening. Note that all characteristics ip,CTp, 7 p 
depend on the wavelength, due to the dependence of the limb-darkening law from the wavelength. 

Till this point, we did not assume that ip is small enough to justify spectra expansions involving powers of ip. Now we 
assume that p is small in comparison with the line widths of ^ then hence ip is small enough to perform such a decomposition. 
In this case it is also legal to process the rotating star spectrum in the way similar to (12.171) . Then we have 


Mo(s) 

MRs) 


= ^(,)_^'(,),p(,)+ 




1 


Tp^(i) +i^(i)] - -.^'"(i)[7p(i) + 3ip(i)cTp^(i) + i^(i)] + ^(p" 


(T2(i) =P 


2 ^2^ 
M*(i) ’ 


M^(i)= j x^I{\R\,s)dR. 
\ R \<1 


(2.19) 


As follows from (12.191) . the stellar rotation does not introduce a systematic Doppler shift or additional asymmetry of line 
profiles in the uneclipsed star spectrum. 

Another method to approximate t hese spectra is to a ssume that they have a simple enough functional shape with some 
parameters to be defined. For example iBoue et al.l (120131 1 use extensively approximations by the Gaussian profile ^p{s — u). 
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6 R.V. Baluev & V.Sh. Shaidulin 


We consider here multiline spectra, so we introduce the following multi-Gaussian function: 

N 

% (■s, «, c) = ^ C 0 p. {s -Ui). 


( 2 . 20 ) 


Whenever deemed appropriate, we may try to approximate 0p or by a function from this family!^ Note that in this way 
of modelling all line profiles become symmetric by definition, whereas (12.181) and (12.1911 may handle asymmetric lines too. If 
we approximate 0 in such a way than from (12.181) we can obtain 


Mo{s) 


1 O-iS'a 

fp{s-Sp{s),u,c) + - (s-Sp{s),U,c)(Jp{s) + . 


N 

/=1 

N 

/=1 


%(■? 


2 '^^(3 

- ip (i) ) + CTp (^) —^ (i - M; - ip (j) ) 


m 


+ ... = 


-(i-M;-ip(i)) + ... (,)(i-ip(i),M,c) + ..., l5li{s)=Pl + a^{s). 


( 2 . 21 ) 


Here we applied an easy identity d'^^{s)lda = /ds^ for each individual line profile. As we can see, formula (12.2111 

reflect nothing more than the line broadening effect by CTp. The approximation (12.211) is valid as far as the multi-Gaussian 
model for 0 is justified, and the decomposition (12.181) is legal. The remaining terms in (12.211) have the same order as in (12.181) . 
namely either <?{xr’) or When neither u nor r is small, i.e. when we deal with a large object eclipsing a fast rotating 

star then the spectrum 0p is not Gaussian even if 0 is. Moreover, its lines might gain signihcant asymmetry and their shift 
might become different from ip. Likely, this case can be only processed numerically, and we do not consider it in our work. 

Based on (12.191) . we may construct a similar Gaussian approximation to the star spectrum: 


M*(. 


y =%(^,'u,c) + ^(u^). 


/3L = A' + ct 2 


( 2 . 22 ) 


but it has more restrictions than (12.211) : it is only legal for small rotation velocities. If this is not fulfilled then is not 
Gau ssian actually, ev en if 0 and ^p are, and to obtain 0^, we should convolve 0 with a specialized rotation kernel, see 
e.g. ( Boue et al. 2ni,lll . The result still might be approximated by a multi-Gaussian function (i, u,c*) with a satisfactory 
accuracy: 

,^*(i)^Mj(i)i#p^(.,u,c*), = + (2.23) 

but there is no guarantee that the broadening parameter CT* is the here same as defined in (I2dll) . although it should be of 
the same order at least. Also, we should introduce the best htting values of line intensities, c^,, which may become somewhat 
different from the original c. 

In our computations we often deal with various convolutions, where the following property might be helpful: 

= (-l)^^^^i^(Mi - M2)- (2.24) 

This identity can be proved by applying a Fourier transform to its left-hand side as to a function of mi — M 2 . In particular, 
scalar product of two multi-Gaussian spectra can be represented as 




/3i 


ds^ 


(j,Ul,Cl)- 




-{s,U2,C2)) = E ciiC2j^‘'''f^^!^{uu-U2j). 

ij=l yPu+Pij 


(2.25) 


In practice we often compare same or close line patterns with ui = 1 x 2 or ui ~ U 2 . In this case (12.2511 can be simplified further. 
If all or the most of spectral lines are well separated from each other (do not overlap), the diagonal terms of (12.251) are 
dominating, while off-diagonal ones can be neglected: 

-(s,U2,C2)} ~ (-!)'" E ciiC2i5^yi?^(Mii - M2,) = (-l)*'ff^^^'"'(0,u',c'). 


d- 

p-(s,Ul,Cl) 


ds^ 




/=1 
_ ft2 


(A0" = i3ii+^2i, 


U = U 2 -U 1 , 


Ci — CiiC2i‘ 


(2.26) 


3 ROSSITER-MCLAUGHLIN ANOMALY FOR A “SMALL” ROTATION VELOCITY 

We consider two types of approximations: 


^ When doing so, we basically subtract the continuum from our spectra. As the continuum is a slowly-varying function, in comparison 
with the lines, its effect on all scalar products like is negligible as long as at least one of A: or m is nonzero. However, the 

continuum becomes important in the norm ||x^||^, so it would be illegal to apply arbitrary normalizations to our spectra without taking 
into account the continuum. 
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Analytic models for the Rossiter-McLaughlin effect 7 

(i) Small transiting planet. This means small planet/star radii ratio and small relative subplanet flux drop f = Mq/Mq. 
However, the rotation velocity t) and hence the subplanet velocity are not necessarily small and may be comparable and 
even exceed the typical width of the spectral lines in Due to the small planet radius, the subplanet spectrum can be 
approximated by (I2T81) . But the expansions (l2T9ll are not applicable. To process this case, we assume multiline spectra 
models with Gaussian line profiles, implying representations (12.2111 and (12.2311 . We consider only first-order approximation in 
/• 

(ii) Small rotation velocity. This means that l) and i-p are smaller than the typical width of the spectrum ^ lines. The 
transiting object (and hence the flux drop /) is not necessarily small and can be comparable in size to the star itself. The 
subplanet and rotating star spectra both can be represented via (12.1811 and (12.1911 . 

In this section we only give our results for the second case, because it is the case in which our results are neat and their 
practical use is easy. The approximation of the first type is considered in Appendix 1X1 


3.1 Cross-correlation with a predefined template 

In the approximation of small rotation velocity we use (l2T9ll to obtain 


^ T 


= i (((Tp2 + i ((7p + 3 ^p(t 2 + ) -f ^(u4), 

i (3.1) 


= 

To transform these expressions to a bit more simple form, we use the property that quantities Mo,ip,Ctp,CT*,yp are all slowly 
varying functions of wavelength, in comparison with the spectra ^ and that contains numerous narrow lines and vary 
quickly. We make an additional assumption that spectral lines are distributed more or less uniformly in the spectral range 
of interest and do not reveal systematic changes of characteristics over the spectrum. In this case variations of any selected 
slowly-varying function A(i) are uncorrelated with variations of ^(s) and justifying approximations of the type 







(3.2) 


This implies that we can just replace the integrals Mi;(s) and M^{s) by their wavelength averages M/, and and define averaged 
quantities ip, dp, d*, yp in exactly the same manner as in (12.1811 and (I2.19II . but replacing Mi; with Mj;. For the sake of simplicity, 
we will omit these averaging overlines from our further notations. Now we can write down: 








(k) 


= Mo 


= Mo 


- .p ) -f ^(ct 2 +.2) _ l(yp + +^3) 


(3.3) 


The wavelength averaging operation on Mi;{s) is equivalent to making the same averaging on the limb-darkening law I{R,s). The 
limb-darkening is usually represented as a linear combination of several simple functional terms that are independent of the 
wavelength, but their coefficients are. Therefore, such an averaging can be reduced to an averaging of only the limb-darkening 
coefficients in I{R,s). 

Now, using formulae &7m and (iTSll . and the constraint (Ireli . we can finally derive an approximation of the RM anomaly: 


= Vi + vV'2 +yVo + > 
fSp 


1 

V = - 




2 ’ 


1 

y = -T 


6 


Vi = - 


V2 = 


l-f 

f 


Vi = 


1 -/ 

/ 

1 -/ 


MiU 

„2 


1-.^ 


M, 


0 -Mo 


M* 2 


Mi 


M*-Mo 


7p + - 


1-/ 


■ +>^11 


1 +/ 


Mi-Mo 


M-i - 3Mi 


+ 


2Mi 


M^-M2 
Mi-Mo ' [M^o-Mof 


(3.4) 


The coefficients V and fl only depend on the star spectrum and on the cross-correlation template. They do not depend on the 
transit geometry and do not vary during the transit. The quantities that vary during the transit are 14■ 

It follows from (iTSl) that in the case when the CCF is fitted by a Gaussian, we may use the same formulae for s as in 
the case of a direct CGF minimization, but replacing by its gaussian-broadened convolution 2^j{s\a). The broadening 

parameter CT should be set to the best fitting value, defined by (Irgl) or approximated in (l2T3ll . Here we must take care of 
the mutual correlation dependence between cr and s. In the formulae (E3, the template is present in the coefficients V 
and fl, so via o they become also dependent on the transit geometry and phase as V(. To solve the task rigorously, we should 
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decompose both a and f into powers of V, solving two equations jointly. Fortunately, such a complicated procedure becomes 
unnecessary. Substituting (ESI in (12.1311 and (I2.1UII it can be easily obtained that the first-order ff(v) term in a actually 
vanishes, so that a = aQ + This means that in the definitions of v and fl we may just replace by ^j{s\(yo), which 

does not depend on the transit geometry again. This would introduce an additional error in drill of which we neglect 

anyway. Therefore, formulae (1321 remain almost the same for the both flavours of the cross-correlation technique. 

To perform a fitting of RV data with the model dSH. we likely need to compute partial derivatives with respect its pa¬ 
rameters, which are necessary for gradient minimization of the chi-square function or other goodness-of-fit statistic. Therefore, 
we simultaneously give expressions for partial derivatives of 14 over and : 


Vi 


V4 


V4 


dVi _ dVx _ MiU 


dV2 

dMo 






1 /■^p 

2MI \ _ 


1 dVi _ 1 _ 1 1 

U dMx ~ ~M^-Mq ~ 1-/’ 


1 1 





2l+/\ 

PI-/; 


SMq dMo dMo Mq ’ dMi SMq ’ dM2 dMi ’ 

\ dV2 _ Mo _ 1 / 

t? dMJ “ (MJ-Mo)MJ “ 1-/’ 

dVo _ dV-i _ / MI-M 2 6M\ \ _ 

mo ~ ~ (M*-MoY M*-Mo ' (M*-Mo)2 ) ~ 


dVo 

JmI 


1 / 


7p + 35 ’j 


(1+/)<^P ”2(72 


1 -/ 


+ ‘^n 


(i-/r 


-3u 


dV2 

dMo 


dM2 dM*’ 


dVj _ n 2 ^^3 _ 2 

dM2 dMo ’ dM-^ dMi 


(3.5) 


The momenta Mj. with their partial derivatives will be computed in the following sections of the paper. 

Although the coefficients V and fj. in (13.41) are expressed by explicit formulae here, we believe that in practice it is difhcult 
to predict them reliably, especially in those works where a reanalysis of public releases of Doppler data is performed, and 
authors do not have access to the full internal characteristics of the Doppler reduction pipeline. In such a case, we suggest to 
treat v and U s-s additional free parameters of the RV curve fit, similarly to e.g. the limb-darkening coefficients. In this case, 
three fittable coefficients of the decomposition should better be defined as 


p' = Pc = Fsin!, v'= VV^c = vvV sini, fl'= sini, 


(3.6) 


because imbed the power factors l)^, and in practice we measure unnormalized Doppler shift Vj. rather than Vj/c. The 
quantities Vt) and are adimensional, while v' and fl' have the dimension of velocity. 

Note that v and fl in (13.41) are dehned via the surface spectrum 3^, which is not observable in practice. Here it is admissible 
to substitute the observable spectrum in place of because the error caused in s by such a substitution is only 
which is neglected anyway. 

Recall that the template should be shifted in such a way as to satisfy (Irel) . which in the approximation of small P 
turnes to 

2 

+ + (3.7) 


We expect that for any reasonably chosen template the value of fl should be positive. This is because we can also rewrite 
6fl = I , and if and have lines in the same or close positions then the both products of the derivatives 

remain positive over the most of the wavelengths range. The value of fl can be negative only if the lines positions in have 
little common with those in or e.g. if there are many emission lines that are wrongly modelled as absorption ones. 

In other words, if the value of fl determined from the observations of the RM effect appeared negative in a particular case, 
this indicates that something is wrong with our model, rendering it unreliable. The coefficient v may have any sign, however. 

Non-zero V may indicate either an imperfect match of the cross-correlation template lines with those in the star spectrum, 
or systematic asymmetry of the line proHles. The latter fact is of a high importance, because this means that asymmetric line 
profiles may require addi tional correction of the RM curve exceedin g, and the order of this correction is larger than of the 
corrections considered by Hirano et al. ( 20ld f and Boue et al. ( 2013ll . The quanitity l/^/Jl is a characteristic of an averaged 
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width of line profiles. For example, for multi-Gaussian and (see App.j^for definitions) we have 

-N -«/// 


V ~ 


y ~ 




3 ^i=i {p^+ii,,y/2 
.(-A.,.) -2 


2i:f=lC,CT,/5 




Aui 


1 fPi+Ph _ 1 W+ftvF" 




)372 


■ + ff{AuA, 


+ e{Avr 


while the constraint (ITfl) can be reduced to 


CiCj.i 


^Aui + Gird) + G{AiA) = 0. 


tx (A'+A?,,)'/' 

Thus, in the approximation of Gaussian profiles, the quanitiy 1/(2^) characterizes an average value for p^ + j5jj. 


(3.8) 


(3.9) 


3.2 Cross-correlation with an out-of-transit stellar spectrum or parametric modelling of the stellar 
spectrum (iodine cell technique) 

Now we should just substitute in place of in the formulae presented above. Formulae (1^ can be transformed to the 
following: 


s = Vi+yV3 + G{v^), 


t^ = 7 


1||^' 


"l|2 


(3.10) 


6 ||.^'||2 ■ 

We can see that now the term with V 2 , which was responsible for either template imperfections or asymmetry of spectral lines, 
disappeared. The coefficient y is now guaranteedly positive. As before, we may treat y as a, fittable parameter of the model, 
if we do not have enough knowledge of the spectra details. We however have a concern that in practice a subtle violation of 
our simplificating assumptions may cause additional disturbing effects in f. Therefore, it might appear reasonable to use in 
practice the full three-term formula (ITil) even for Doppler data obtained with the iodine cell technique. At least, it might be 
a matter of practical verification with real data, whether the term with V 2 indeed becomes negligible in this case. 

For multi-Gaussian ^ we obtain 


A I- ft's 


N .2 


In this approximation, 1/(4^) measures an average value for jS?. 


(3.11) 


4 COMPUTING THE RV MOMENTA 

We need to compute the integral momenta 

Mk[d,r,X)= I Al{\R\)dR, M^=Mk\s^o^r=l= / Al{\R\)dR fc = 0,1,2,3. 


(4.1) 


S{S,rA) |ii|<l 

Here we consider them as functions of three parameters 5,r,2,. Their definition and geometrical layout are given in Fig. [T] 

The star radius is assumed unit. _ _ 

To compute these momenta, we apply an adaptation of the approach developed by Abubekerov fc Gostev (l2013h for 
transit lightcurve modeling. Dehne auxiliary functions 


y/x, JC > 0, djuf 

0, X < 0, dx 


= -^ 


l-x2 


'V{d,x,y)=x/ 


/ 52 -|-x2 —y 2 
y 2x5 


i n, x<-l, 

arccosx, |x| < 1, ,S(x) ='3i^/x = 

0, X > 1, 

(4.2) 

Then, extending formula (7) from ( Abubekerov fc Gostev 2013h to take into account additional factor x^, we can write down: 
1 'V{S,p,r) 

Mk= I{p)p'‘^^dp j cos* ((P +A) dtp. (4.3) 

b -'!'(6,p,r) 


Abubekerov fc GostevI (l2013l ~) computed Mq and its derivatives with respect to r and 5, which are necessary to express the 


flux reduction AL and its parametric gradient (for further htting of the model by a gradient descent). Here we also compute 
M]^ with their derivatives for k= 1,2,3. 
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Figure 1. Geometry of the transit illustrating the Rossiter-McLaughling effect. 


Let us put p = sin6. The limb-darkening law is often modelled by a polynomial in cos 6. Then can be expressed via 
linear combinations of the following integrals 

f A+'P 

Hnk = y^sin*^* 6cos”^* J cos^ (pdip, T ='i'(5,sin6,r), 

0 A-'P 


Assuming at first that k is odd, rewrite cos^ (p as a trigonometric sum of multiple argument and compute the inner integral: 


2 A+'P [k/l] 

y sin^^* 6cos”^* ^ cos(fe —2y)(pd(j!) = 


I I [k/'A 

= —j^ / sin*^+* Scos”^'6 ^ ^ sin(^ —2y)*i'cos(fc —2y)A = 

I [k/A } 

= ^ —^cos(fc —2_/)A / sin^^'6 cos”^'6f4_2j-l(cosT)sinTdO. (4-4) 

2 j^ok-2j J 


where (]„ are Chebyshev polynomials of the second kind. If k is even, then the procedure is similar, but the sum should also 
contain a term with j = k/2, which must be handled separately due to a degeneracy. In general, we can represent H„k via the 
following trigonometric polynomial in X: 

[k/A 

Hnk= ■^nkiCos{k-2j)X, (4.5) 

./=() 


where 

£ 

j 2 

’^nkj = ysin''+'ecos"+'ef/i_2f-l(cos'T)sin'Tde, 2j<k, 

^ 0 




hi- 


6i 




e'Vde. 


(4.6) 
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By extracting from cos"+* 6 a multiplier cos^ 6 = 1 — sin^ 6, we can obtain the following recurrent relation: 

■^nkj = , 2; < k. (4.7) 

This can be used to reduce the index n to either « = 0 or « = 1, by the cost of increasing k and j. In practice we only need this 
formula to reduce n = 2 (quadratic limb-darkening term) to n = 0. 

Now let us define 


fl=l —(5 —r)^, a(6) = a — cos^ 6, 
This yields 


cos Tsine = (5 — r) + 


a{6) 


sin'T sin 6 = 


6 = (5 + r)^ —1, b(6) =b-\-coA 0. 

^[a{e)b{e)] 


(4.8) 


(4.9) 


28 ' ... 28 

In fact, the integration range can be limited to only the range where sinT ^ 0, corresponding to —b < cos^ 6 < a < 1. For 
convenience, we need to transform this variable range to a constant one. The case a < 0 is trivial, as sinT = 0 everywhere. In 
the case a> 0 we can introduce the following replacement 6 i-> f: 

cos^ 6 = a — min(a,a-|-6)r^, sinO cosO 86 = min(a,a + b)t d/. (4.10) 

After this replacement, the new integration variable t always spans the same segment [0,1]. By making this replacement, 
expanding Chebyshev polynomials as I/«(v) = ^-nd introducing additional auxiliary designations, we may rewrite 

the integrals .J^nkj as follows: 

1 




nkj 


Ck-jK 


cj p 

2^-2 k-2j28 

u 

k-\ 

x{q^+x^ ■ Uk-2j-\ 


JPk-j {pd) Y^max(l,m) — f2 ^1 — min(ra, l)t2^ ^ dt, 2j<k, 

q+J5 




rtii 


= v ^ Uk-2j-l,l 
1=0 




X \k-2j-2l-l 

2S) 


{q^+x) 


j+l 


q = 8-r, 


Ar8 


p = min(4r5,1 — ^ ) = 4r5min ( 1, — ) = (1 — ^ )min(l,m). 


(4.11) 


Note that although there is a division by 8 in p/{28) and x/{28), in actuality there is no pecularity at 5 = 0, because by 
definition p < 4r5 and thus p/(28) < 2r. Note that 0 < m < 1 corresponds to full phase of a transit (r + 5 < 1), while m > 1 
corresponds to a partial occultation |r—5|<l<r+5. The cases of a full eclipse or no eclipse (|r — 5| > 1) would correspond 
to m < 0, which are illegal in these formulae by definition. In the latter case, ^nkj = 0 for 2j <k. 

The case k = 2j is more complicated due to a “naked” T in the integrad. We may apply two ways of integration by parts: 




5cos"+2e „ 2n&(-q) 2/-1 ^ 1 C; 

“ = ^2 „ + 2 + -i^^n+2,2j-2J-l + 




2 

A 


'¥cos"e 


+ 2 

5sin22+2e 

2j + 2 


n + 2 


1 + 2 221- 


,/sn 


— ; ■ . 1 


n 1 C; 

^^■A-2,2J+2J+1 j 


1 C2J f . 
1 + 12217 ®* 


A 6005"+^ e^de = 
dO 




sitdj+^ecos"e—de, ( 4 . 12 ) 

du 


221-17 2 ;- + 2 ^"221 ; + l 

0 " ■ 0 

where 5/j; is Kronecker delta (not to be mixed with the planet-star distance 5, which is unindexed), and we use that 
limp^+()'i’(5,p,r) = 7l&{—q) (with 0 being the Heaviside function). Taking into account (14.71) . both formulae (14.1211 appear to 
yield the same recurrent relation, if n ^ 0: 


[n + 2j+ 2 )^fj 2 jj — (2j A'^n, 2 j— 2 q—\A-Gfi-\- 2 j-, J > 0, 


0 


^2 < sird-j e cos" 6 A-dO. 


dO 


(4.13) 


For n = 0 this relation is still valid thanks to the first of eq. (I4d^ . but we also obtain from the second eq. of Km an 
independent non-recurrent formula 




'I'(5,l,r) Gqj+i 


j>0. 


(4.14) 


0.271 221 j+l 2j+l’ 

It can be verified by direct substitution that this formula satishes (14.1311 . The recursion (14.1311 can be used to reduce j until 
we meet j = 0. For j = 0 we may use one of the following schemes: 

(« + 2),y„,O.o=2;r0(-?) + G„+2,o, ^o.Ofi = 'V{8,l,r)-Go^i. (4.15) 

We use the second formula of (14.151) for n = 0, because this simplifies some computations. 
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To compute Gnj, we express the derivative of T by differentiating cosT and then apply the substitution (14.1011 : 


Gnj = r[\- 




2\>-l (l — min(l,OT)r^)- 


max(l,OT) 


- dt. 


(4.16) 


Note that pl{2r) < 25, so there is no pecularity at r = 0. Note that in the case of a full eclipse or no eclipse, |r-- 5| > 1, we 
have m < 0 and put G„j = 0 by definition. 

Now let us review the results. Our final useable formulae are (14.514.714. 11 14.1314.1514.1^11 . For even n all integrals are 
elementary, because their integrads are rational functions of t and of the radical . For odd n, the integrands 

are rational functions of 1 and (1 —min(ra, implying that all integrals are elliptic and can be expressed via 

Legendre complete elliptic integrals with the parameter min(m, 1/m). In the most cases, this should be only the Legendre 
integrals of the first and second kind, and only in Gi.o we meet the Legendre integral of the third kind, where t is present in 
a denominator of the integrand’s rational part. This integral, however, affects /i due to the recursion (14.1311 . 

Our approach allows to give exact and explicit formulae for all velocity momenta for any polynomial limb-darkening 
law. However, in this work we only need momenta up to cubic order {k = 0,1,2,3) and up to quadratic limb-darkening law 
(n = 0,1,2). From now on, we stop using generic notations and focus on the computation of Hnk for the specified indices n and 
k. 

Consider the quadratic limb-darkening model 

/(p) = + (A, 4-2A,)(/^ -/) +A,(2/‘ -/^), 


/‘=l-COSe, /' 

Therefore, using (iTTl) for terms with n = 2, we have: 

Mk = (1 - A, -2A„)M| 4- {Ai+2Ag){M^^-Ml)+Ag{2Ml-M^^) 

Mq = = "A ,0,0, 

M^ = JPiqcosX, Mf—M\ = J^iiocosA, 

= J^o,2,OCOs21 + 3^0,2,1, “^1,2,0^08 21 + ,.^1,2,1, 

= J^o, 3 ,ooos 3/1 -I- 3^0,3 ,1 cosl, = ,.^130 cos 31 + ,.^^ 3,1 cosl 


T = cose, /^ = (i-cose)% = i-cos^e. 


1 = 0 , 1 , 2 , 3 , 

2M> - m;] = 2J?b, 2 , 1 , 

2MJ-M5‘ = |j^O, 3,1 cosl, 

2 M 2 ^1^2 = “.^0,4,1 cos 21 + ^,.^0,4,2, 

2 M 2 —M 2 = 5“^0,5,1 cos3l + I J^o, 5 , 2 COsl. 


(4.17) 


(4.18) 


Thus, we need to compute 16 integrals in total: 11 of them are of type k> 2j (with 7 elementary, and 4 elliptic) and 5 
are of type k = 2j (3 elementary, 2 elliptic). Note that 

-^0,0,0+ ^2,1 ^ _ 3J^o,2,i+G 2,2 _2;r0(-^)+ G3 ^o 


Ao ,0 = 'I'(5,l,r)-Go,i, 


0 , 2,1 = 


0,4,2 = 


“1*^1,0,0 = 


“1^1,2,1 = 


“1^1,0,0 + <J3,1 


(4.19) 

The remaining part is to compute integrals (14.1114.1611 for the indices specified above. This is a routine but difficult work 
due to quickly growing formulae. We used MAPLE computer algebra to compute these integrals in a symbolic form. The 
results are given in Table [T] We represent all y„kj and their derivatives as linear combinations of the following functions: 
elementary y/,(p, and Q for n = 0,2 and elliptic A,E, and H (or f) for n = 1. The coefficients are algebraic polynomials in r and 
S. The degree of these polynomials may reach 8, and we tried to reduce them by reusing an auxiliary function W, whenever 
possible. All these analytic expressions were verified by comparison with numeric calculation of the integrals. 

In this work we adopt the following definition of the complete elliptic integrals (Legendre normal forms): 


K{m)= [ — 

i ^ 


dO 


- OTsin^ 6 


E{m) = j V 1 — rasin^ Odd, n()i|ra) = J — 


dO 


• 6 )\/1 —OTsin^ 6 


(4.20) 


0 '' 0 _ 0 __ 

Note that the sign of n he re is opp o site to what is adopted by ( Abubekerov fc Gostev 2ni3h and by MAPLE, but it coincides 
with what is adopted by Carlson ( 1994h and by the GNU Scientific Library. This choice allows to keep th e argument n 
always positive in our formulae. To compute these integrals we recommend the algorithms by Fukushima ( 2013I L This method 
appears faster by the factor o f a few in comparison with the lCarlsonI (1199411 s ymmetric forms a pproach, which was selected by 
Abubekerov fc Gostev ( 2013h and also adopted in the GNU Scientific Library. Fukushima ( 2013ll uses the following “associated” 


forms of elliptic integrals (note that we also changed sign of n here): 

EE E 

cos^ 6 d6 , / sirP'6 d6 


, , f cos 6 da , , /■ sii 

SH = / , . 2 , D{m) = / — 

/ V 1 — sin^ 6 n V 1 


OTsin^ 6 


J(n\m) = / — 

■I ( 1 - 


sirr Odd 


0 '' O'" 0 

This implies: 

Kim) = B(m) A D(m), E(m) = B{m)A(i—m)D{m), Yl{n\m) = K{m) ~ nj(n\m). 
Additionally, the following identity from ( Fukushima 2013ll becomes useful for us: 


0)2/1 — msivd 6 


J{n\m) = -- 


2 ^n{l+n){n + m] 


1 . , m fm\\ 

H— Kim) [ — \m) . 

n \n\ / 


(4.21) 


(4.22) 


(4.23) 
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This identity was used to remove an undesired discontinuity near 8 = r that appeared in the original MAPLE output, which 
was expressed via n (see intermediary quantity Q. in Table [T]). 

Some formulae in Table [T] contain an apparent pecularity near 8 = 0 due to division by 8. All functions are actually 
smooth at 5 = 0, but the pecularity is associated with the subtraction of close number like KR) — ER) for e ss 0. This leads 
to an accuracy loss near 5 = 0. To get rid of it, we might rewrite the formulae using some non-standard elliptic functions as 
a basis, but this is not convenient, so we choose to consider the case 5 ss 0 separately and provide in Table [2] the relevant 
Taylor series about 5 = 0. These series are more preferred than general formulae, if 5 < 0.05(1 — A). We give these series only 
for elliptic integrals n = I and only for the case of the full phase of a transit, r + 5 < 1. The case of a partial occultation is 
possible with small 8 only if an additional condition r ~ 1 is satisfied. In practice it is a rare condition when simultaneously 
8 is small and r is close to unit (a close-to-ring eclipse). Besides, numerical tests did not reveal significant loss of precision in 
this parametric domain. 

_ Note that our results for Mq (which is responsible only for the flux decrease) are in agreement with lAbubekerov fc Gostev 

1 2 OI 3 II . Also, substituting 5 = 0 and r = 1 in Mj, we obtain values for the whole-disk momenta 


M^ = k- - A, - -A. 


TZlTZ n 

MX = -A;-Ao, 

2 4 60 15 ^ 


3 = 0. 


(4.24) 


3“' 6'^’ 

In addition to the expressions (l4T8ll we may also need to compute partial derivatives of with respect to A, which are 
trivial, or with respect to the planet coordinates x,y in the projection plane. We can use the following formulae for this goal: 


r.9 9 3 88 A uu y 

8^=R+/, — = -, — = 4 

ax 


X 

8 ’ 


dy 


y_ 

5’ 


ink = T„ 




5 cosuA 
dx 


= n^U„^ 


ii)- 


d cos nk 
dy 


^ rr 


(J)- 


(4.25) 


Here and f/„ are Chebyshev polynomials of the Erst and second kind. 
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Table 1. Integrals ^„iij{S,r) and their derivatives: general formulae. 


►1^ 


Function / 

'^000 

^010 

.f^020 

d^030 


.1^021 

d^031 

■Ami 




J%42 


•Aoo 
•A 10 

-^120 


-^130 


.-^121 


-<^131 


Auxiliary 


/(5,r) 

. '> o 

\lf + r<p-^ 

r25<p + (l-52-r2)^ 
r252f + [352(l-r2-52)-W]4r 
1+ [652(/ - 3r^S^ - 2r^ -S^ + l)- 

^ + r^{r'^ + 25^) f - (5r^ + 5^ + 1) 

r25(r2 + 52)^ + [3(r2 + 52)(l-r2-52)-21F] 

r252(3r2+ 25^) I - [(r^-35^ + 3)133+ 

+6S^r‘^ + 9r^S^ + 25^ + r^- ^2)] ^ 

(2r^ + ^ + [305“ (2r2 + S'^)(\-r^-S'^y 

-(l + 55“ + 252-r2 + r252)5W-lF2j 2_ 

|+r2(/ + 6r252 + 35“)|- 

- (10/ + 5“ + 19/5^4/ + 5^ + 1) ^ 

'48 

/5(/ + 3/5^ + 5“) ^ + [(3 + / + 8'^)W- 
8 

-6(5/5“ + 5/52 + 3/52 + / + 52-1)] ^ 
|n+| [(7/ + 52-4)£-IF^] 

ifj [(16/52-lF)£ + (l-52-/)vi3^] 

{ - [5^(9/ + 75^ - 7) + 21F] WK+ 

+ [8/52(/ + 155^ - 1) + (2/ - 55^ - 2)1F] E} 

{- [52(224/5^ + 6352 - 64/ + 127/ - 63)+ 
+(8 - 8/ - 3552)1+] WK+ [l6/5“(8/ + 7252 - 9)+ 
+ 52 ( 29 / - 2752 - 36)1+ + 81+2] E) 
^i2-^[(39/+952 + l)l+^- 

- (1 29 / + 95“-68r2+ 246^2 52 _ 3^2 _ 3 j 

- ^ {( 32/52 + 7 / + 752 _ 2 _ ■iw)^/^K+ 

+ [16/52(1 - 8 / - 852) ^ (4_|_ 3^2 3 _ 352 )^j g| 


dfjdr 

2 r(p 

2 rS ( p - ^ 

r52(p + (/-352-1)^ 

r521 + [21+ + 352 (/ - 352 - 3)] 

r(/ + S ^)( p-rQ 

r5(2/ + 52)^-(l + / + 552)|« 
r52(3/ + 52)9+ 

+ (/ - 175“ - 8/52 + / - 552 - 2) ^ 

r52(4/ + 52) ^ + [(3 + / + 3752)1++ 

8 

+652(7/-17/52-13/-1452+4)] 
r(/ + 4/52 + 5“) ^ - (1 + 3/ + 352) ^ 

r5(3/ + 6/52 + 5“)^ + 

4 

+ [1+ - 3(7/52 + 35“ + / + 252)] ^ 

4r£ 

If [(/+752- l)£-l+^] 

^{(/-452-1)1+^+ 

+ [ 52 ( 4 /+ 2052 -5)+!+]£} 

{ [352(5/ - 2152 - 14) + 81+] WK - 
- [352(104/ - 232/52 - 17552 - 209/ + 105)+ 
+ (8/ + 28152-8)1+]t} 

If [(4/ + 452- l)F-l+^] 

A [(88/52 3_4Qg4 ^5^2 _ 552 _ 5 _ j,w)E- 
lOo '■ 

-(3/+ 1352 + 2)1+^] 


df/dS 

_Q 

/9 + (/-352-l)A 
/5(P+[1+-352(/ + 52 + 1)] ^ 

/52 ^+ [252(3/ - 9/52 - 4/ - 552 + 1)+ 

+(1 -/+ 552)1+]^ 

/5(p-(/ + 52 + l)^ 

/(/ + 352) ^ + (/ - 55“ - 20/52 + / - 552 - 2 ) j ^ 

/5(3/ + 452) I + [(/ + 952 + 3)1++ 

+652(/- 15/52-5/-452)]^ 

/52(6/ + 552) ^ - [(/52 - 95“ + 3/ - 2452 - 3)51+ - 31+2+ 
16 '■ 

+3052(/52 + 21/5“ - 5/ + 14/52 + 8/ + 752 - 3)] 

/5(/+ 52)^+ [l+-3(4/52 + /+ 52)] g 


/(/ + 9/52+55“)^- [(/-752 + 3)1++ 

8 

+6(7/52 3^ 23/5“+5/52+45“ - / - 52+1)] 

^[(/ + 52-1)£-1+^] 

{(/ - 452 - 1)1+^+ [52(19/ + 552 - 5) +1+] E } 

{[41+ - 52(17/ + 2152 + 14)] WK ~ [(4/ + 2552 - 4)1++ 

+52(16/ - 3552 - 320/52 - 51/ + 35)] Ej 

( [5“ (40/ + 1368r2 52 + 39^2 _ 105 _|_ 10552) 

630d'^ ^ '■ 

+52(13/ - 7852 + 36)1+ - 81+2] £ + [(8 _ g/ + 7052)1++ 

+52(83/ - 259/52 - 125/ - 14752+42)] 1+^} 

jIj [(48/52 + 5/ + 552 - 5 - 31+)£ - (3/ + 352 + 2)1+^] 

{ -(59/52 + 215“ - 7/ + 752 + 7 + 31+)!+^+ 

+ [52(152/ + 488/52 - 19/ + 3552 - 35) + (4 + 3/ - 1852)1+] E } 


i// = T(5,l,r)€[0,;r], (p = 'I'(5,r, 1)€[0, tt], a = l — (5 — r)2, 5 = (5 + r)2 — 1, W = ah, Q = ^[1+], Q nonzero only if r + 5 > 1 and |r — 5| < 1; 

^ = £ = -K ifr + 5 < l (implies |r-5[ < l); 

K = K {^)/ V ^, E = V ^ E (^)- bK , n = f ^( 4 rS )- y ^ j[^,^)-K if r + 5 > 1 and |r-5| < 1; 

^ = £ = n = 0 if [r-5| > 1 (implies r +5 > 1); f2(5,r) = ;r©(r-5) + |4|(n+(/- 52)^ = J + (52 -/)/ (see text for proof). 


R.V. Baluev & V.Sh. Shaidulin 



















MNRAS 000, [TH23] (2015) 


Table 2. Integrals and their derivatives near 5 = 0. 


function / f(S,r) 


3f 


3f 


-<^100 

¥[i-< 

.1^110 

TzdSVl 


-^130 


Am 


2KrVY^ [l + ^ f 1^] + ^(5^) 


2 nrS 


VT- 




3/-8r2 + 8 52 
2 (1 —16 
r'>-6/ + 24r2 + 16 S* 


rtT, [1 5 „2 


nrS^ 

(1_ ^2)3/2 


128 


+ - 


(1_.2^4 
!r- 

5/- 18/ + 24r2- 16 5^ 


+ ^(5^) 


l-SZ + ^rV 


nr^S^ 

4(1 — r^)2/2 


( 1 - 


32 
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+ ^(5’) 


37 irS 
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3 4 '' (l-/)2 32 
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45 / - 114/ + 88r 


+ ,^(5S) 

*-16 52 


(1 - /)'* 


64 


+ ^(S") 


(1 — r 
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256 




nY\pi- 


+ ^(5’) 


1 + 


3/-4 352 
(l-/) 2 ^ 
/-4/ + 8 55'* 



TCrS^ 

r, 9 2 45 4 35 6 

3kYS^ 

[^-2’-" + y/+ 

2(1-/)V2 

1 —/h-/-/- 

2 8 16 

4(1 — /•2)V2 


(l-/)4 64j+^(^^) 

nP-b [1 5„2 5r‘^-12r^+8 6^1 , 


5 2 4 

+24'' + 

35 / - 120 / + 144/ - 64 552 


(l-/)2 384 

11)^5 (i 3 ..2 9d-22++16 52l , 

2*^ (l_+)2 16J+'2^0J 


+ ^(5«) 


3;r 


■ \/1 — / / + 


15/-24/+ 8 352 


(l-/)2 8 

15/-52/+ 64/-32 55“ 


(1 - / 


64 
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Assuming that r + 5 < 1 (full phase of a transit) that implies r,5 < 1 and [r—5| < 1. 


Analytic models for the Rossiter-McLaughlin effect 15 













































































16 R.V. Baluev & V.Sh. Shaidulin 



log(wavelength) [arbitrary units] 


Figure 2. Synthetic star surface spectrum used in the test simulation (Sect.[5]l. The spectrum contains 100 Gaussian lines with random 
characteristics. 


5 INVESTIGATING MODEL APPLICABILITY BY NUMERICAL TESTS 

As far as our model employs power-series decompositions in u, it should be applicable only if p is below some limit. How¬ 
ever, such a limit is difficult to predict theoretically. To determine the actual domain of the applicability, we use numerical 
simulations. 

First of all, we construct a simulated spectrum ^ containing only Gaussian lines with randomly chosen characteristics 
(Fig. [21). This spectrum we use to perform direct numerical integrations in (12.11) . We compute the out-of-transit spectrum 
on a grid of V, and the in-transit ones on a 4-dimensional grid of the parameters V,r,S,X. After that, we fit numerically 
each in-transit spectrum —with a shifted determining the best fitting shift. In such a way, we obtain a set 

of simulated Doppler shifts as a table function of the gridded values (p,r, 5,A). Simultaneously, we compute our analytic RM 
model Isa on the same grid, and then compare the results. 

However, it is very difficult to seize a four-dimensional space, so we need to convolve some of the dimensions somehow. We 
consider p as our primary parameter of interest, and for each selected P we compute only r.m.s. of the RM model residuals, 
corresponding to (r, 5, A) from the grid. The domain of the grid was constructed as follows: r e [0.05,1.22], S G [max(r— 1,0), 1 -fr], 
A G [0,;r]. The values of P were ranged approximately from 1/3 to 10 times the average line width. A subtle but important 
detail in this algorithm is that different r imply different amplitudes of the RM curve, scaling roughly as R. Therefore, we 
introduced an additional descaling factor of 1/r^ to equibalance the contribution of different curves in the cumulative r.m.s. 

To compute the RM model we apply three distinct methods to determine the coefficients v and fl- In th® first case, 

we derive them from the original spectrum which in practice would not be accessible to the observer. In the second case, 
we derive v and fl from the observable spectrum And finally, in the third method we just fit our simulated RM curves 
with (IXH) . assuming that p, v, and fl are all free regression coefficients. As we discussed above, the first two methods are 
equivalent for small p, whenever only the three decomposition terms are significant. But for larger p more terms enter 
in the game, making it important, which of the spectra was in use for v and fl- The third method is the one that can be most 
easily implemented in practice, if the original spectra are not available at all. 

Our main results of the simulation are demonstrated in Fig. [3] In this figure, we have four plots that show how the 
following quantities depend on P: (i) normalized r.m.s., (ii) relative bias of the best fitting p, (iii) predicted and best fitting 
values of v (all must be negligible in our model), (iv) predicted and best fitting values of fl. From these plots, we can draw 
the following main conclusions: 

(i) It is much better to use for computation of the coefficients v and fl, while using infers larger errors. It is very 
favourable to us, because the spectrum that we obtain from observations is also and not 

(ii) If we use to compute v and fl, the maximum value of the rotation velocity, after which our model becomes inaccept- 
able, is only about the average linewidth (for the lines in the original spectrum If we use to compute V and fl, this 
limit increases to 2 — 3 times the average linewidth. 

(iii) The best fitting values of v and fl are close to those derived from 

In the summary we may note, that the range of applicability for our model appears rather optimistic. Even though we 
used spectra decompositions into powers of p, our model remains accurate even if p exceeds the average linewidth. In fact, it 
is more adequate to say that our model requires that “p is not so large” instead of “p is small”. 
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Overall model quality 


Relative error of the best fitting u 



Coefficient v (expected negligible) 
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1 

u, relative to average linewidth 



Coefficient jj. 


10 



Figure 3. Results of the test simulation from Sect. The abscissa is the value of t) normalized by the average linewidth for lines in 
the original spectrum ^ (not the rotation-broadened The ordinate in the first plot is a normalized value (r.m.s.)/u. The value of 
r.m.s. also contains an internal normalization by l/r^ (see text), so the graph virtually shows an average relative error of the model (i.e., 
residuals are normalized by the amplitudes of the relevant RM curves). In the last two plots we show adimensional quantities Vl) and 
tlt)^, which reflect the relative contribution of the corresponding correction terms in the model II3.4II . See text for details and discussion. 


6 PRACTICAL APPLICATION: THE CASE OF HD 189T33 

As the number of formulae appearing above was large, let us now describe a concise step-by-step scheme to compute the RM 
anomaly: 

(1) Compute J^nkj based on formulae from Table[T]or[21 if 5 < 0.05(1 — R)- Depending on the expected degree of RM anomaly 
approximation (1,2, or 3) and degree of the limb-darkening model (0,1, or 2), not all of these 16 integrals may be actually 
needed. Whenever the RM anomaly curve is not plainly modelled but also fitted based on the RV data, also compute partial 
derivatives of J^nkj to be used in the gradient minimization of the chi-square or other goodness-of-fit function. 

(ii) Compute momenta from (14.1811 . Again, not all of them may be necessary, depending on the particular practical task. 
If needed for further fitting of the anomaly, simultaneously compute the derivatives of M]^ with respect to x,y, and r, based on 
eq. Km . 

(iii) Use formulae (El) or (13.1011 to finally compute the RM anomaly. If necessary, the gradient of the model with respect 
to the parameters can be obtained based on (13.51) and on partial derivatives from the previous steps. The coefficients near 
Ui, 2.3 can all be treated as free parameters of the fit. 

During the transit, the quantities d and A are varying along the projected planet trajectory. This algorithm is not 
responsible for modelling the planetary orbital motion during the transit, which must be carried out separately, e.g. based on 
a Keplerian or A-body model. We omit a consideration of such models in our paper, as this topic is already investigated quite 
well. 
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Figure 4. Light curves (1 —/) for several sample transit configurations shown on top, and plots of the RM anomaly terms Vi, V 2 , and 
V 3 exposed in successive rows downwards. In each graph, the abscissa is a normalized time with first contact at —1 and fourth contact 
at +1. The scale and labels of the ordinates are omitted, because we intend to demonstrate only shapes of the curves here. The plots 
assume the limb darkening with A; = A,, = 0.25. 


We do not take into account the effect of finite light speed that may cause a subtle time delay between the RV variation 
due to the planetary orbitaly motion and the RM anomaly. This delay appears because the former RV shift is imprinted when 
the light is emitted from the stellar surface, while the latter one appears when the light is blocked but the transiting planet, 
which occurs closer to the observer. This delay should be usually small, e.g. ~ 10 — 20 sec for a typical hot Jupiter. This effect 
is not very hard to model, but this falls out of the scope of the present paper, so we neglect it. 

First of all, let us gain some impression of the behaviour of the basis functions ^ 1 , 2 . 3 - During the transit, they can be 
viewed as functions of the time, and they also depend on the transit geometry. We plot them in Fig. [4] for several sample 
cases of the planetary orbit orientation. As we can see, the case in which the planet motion is parallel to the projected star 
rotation axis, is degenerate. In this case all ^ 1 , 2,3 have the same or almost the same shape, so it would be impossible to fit the 
relevant coefficients separately. But in other cases the shapes of the basis functions are different, and their coefficients can be 
fitted independently. 

Now we apply our RM effect models to the transiting planet of HD 189733. This is indended to be just a preliminary and 
demonstrative study. Full analysis of this and other objects with the RM effect is prepared for a separate work. HD 189733 
is studied very well already, and it offers an ideally suited a testcase. We use T ERRA jAnglada-Escude fc Butl^ l2012l 'l 
Doppler data derived from the HARPS and SOPH IE spectroscopy and published in I Baluev et al.ll2015ll. Additionally, we use 
public Keck RV data given by Winn et al. ( 2006h . and public transit photometry from jBakos et alj bood : Winn et al. 2007 : 


Pont et al.ll2007l '). We do not use a few HARPS-N measurements of this star from llBaluev et al.ll2oi5fi . because they appeared 
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entirely erratic after a closer look (this is being i nvestigated). Also, we do not use vast photometry available for this object in 
the Exoplanet Transit Database, as was used in (IBaluev et al.ll2015l ). Including this photometry slows the computations down 
dramatically without making significant changes to the models of the RM effect in Doppler data. 

We split HARPS data in 5 independent subsets, corresponding to 4 transit series and 1 out-of-transit one. The Keck 
data were split in two subsets, corresponding to 1 in-transit series and to the remaining randomly distributed measurements. 
Finally, three Keck points that were obtained before its CCD upgrade in 2004 were removed. The splitting in such subsets is 
necessary because the RV scatter on the timescale of a single in-transit run is only about ±2 m/s, but on larger transit-to- 
transit timescales it in creases to 10 — 20 m/s. This is likely an activity-induced red noise effect similar to the one considered 
e.g. by ( Balueviboidbl) . In the our case of HD 1 89733 it is easie r to assign fittable RV offsets to different in-transit runs instead 
of dealing with correlated noise models as in I Baluev 2013bb:_^l Dopp ler and tra nsit data were tra nsformed to the same 
BJDxdb time system using the public software by Eastman et al. ( 201 fll) . Note that Winn et al. ( 2006ll published their K eck 


fealuev et al. 201!Th 


data without performing the barycentric reduction having amplitude of ~ 4 min, and RV data from 
in the UTC system, which currently differs from TDB by approximately 1 min. Such differences of a few minutes become 
important for self-consistent RV-ftransits fits. 

In Fig. [S] we show RV residuals corresponding to different models, in the vicinity of the transit, and separately for the 
HARPS, SOPHIE, and Keck datasets. We can see that the RM effect is obvious, its curve is well sampled and measured with 
high accuracy both by HARPS and Keck, whereas SOPHIE has only sporadic data in this phase range (top raw of plots). 

In the second raw we plot RV residuals of the classi c RM model with th e limb-darkening c oefficients fixed at A * = 0.65 and 
A*^ = 0.15. These values are close to those adopted by Triaud et al. ( 20091) . And confirming I Triaud et al. 2001)1) . the classic 
model of the RM effect leaves certain systematic wave-like variation in the residuals, which is clear in HARPS and less clear 
but noticeable in Keck data. Note that all 4 HARPS transit runs are plotted over each other in a single graph. Their systematic 
variation cannot be due to the effects like asteroseismologic oscillations, which would change the shape from one transit to 
another. This variation definitely reflects an inaccuracy of the classic RM model. 

Formally, this variation can be equally fitted by either (i) adjusting the limb-darkening coefficients or by (ii) using the 
correction terms of dsa. These ways offer practically equivalent models. The residuals look almost identical for these fits, and 
in the both cases they leave no significant hints of any other systematic variation (third and fourth raws in Fig. [5]). However, 
in the case (i), their estimations of the limb darkening coefficients appear too different from the theoretically predicted values, 
and actually do not look physically reasonable. This indicates that the case (ii) is more realistic. In this case we obtain a 
well-constrained estimations of the coefficients v' and fl'. The value of v' is close to zero, while fl' appears comparable to V sini 
(see labels in Fig. IS). From (ETH), this value of y’ corresponds to the average width of the spectral lines of ~ 2/3 of V sini, or 
~ 1800 m/s (before the rotational broadening). 

The estimation of V sini in the case (ii) is reduced by ~ 25 per cent in compar ison with the case (i). This reduced value 

with the classic RM model (~ 3300 m/s) 
based on the line-profile tomography. We 
need to emphasize that our model is sensitive t o the adopted v a lues o f the limb-darkening c oefficients, and by incr easing of A* 
we would obtain an estimate of V sini closer to Cameron et al. ( 201ol ). From the other side, Cameron et al. 1 2010|) use only a 
linear term in their limb-darkening model, so at the current stage it is still unclear, which of the two latter estimates is closer 
to the truth. It is however definite that all values of V sini that rely on the classic RM model are overestimated. 

We also tried to fit simultaneously the RM correction and the limb-darkening coefficients (fifth raw in Fig. [Sjl. In this 
case we obtained an ill-conditioned fit with large uncertainties, and the residuals did not change. However we point out that 
the coefficient v is always determined robustly and with a good accuracy, so it does not seem that there is large correlation 
between v and fl or between v and the limb-darkening coefficients. Also v is always consistent with zero within narrow limits. 
This is exactly what theory predicts, because all Doppler data that we used here are obtained by means of the spectrum 
modelling rather than by correlating with a template. In fact, a zero value of v indicates that these spectral models are of a 
perfect quality. 


of ~ 2900 m/s is significantly smaller than the one obtained bijTYiaud^_et_^lJ ([200 
and even smaller than the value of ~ 3100 m/s, obtained bv ICameron et al.l ( 201 


7 CONCLUSIONS AND DISCUSSION 


This paper represents an attempt to construct more general but still analytic models of the RM effect with a particular focus 
to an improved practical usability, especially by a third-party analysis work. Although our primary new model (13.41) does not 
depend on several important restrictions, like the single-line spectrum, or specific line profiles, or small planet, there is still 
much to be done in this topic. The main vulnerability of this model is that it relies on decompositions in V sini, r equiring it to 


be sm all. In fact, we considered both modelling approaches: employ power-series decom positions in V sin i, as in (IHirano et al. 


2O10l) . or avoid such decompositions by assuming a simple Gaussian line profile, as in ( Boue et al. 20131 ). However, our most 


useful results correspond to only the first case. In the seco nd case we did not succeed very much, just showing that the actuality 
is significantly more complicated than explicated e.g. by Boue et al. ( 20131 ). Nevertheless, we believe that our primary model 
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HARPS SOPHIE KECK 

Only modelling the orbital variation 



orbital phase orbital phase orbital phase 


(orbital variation + classic RM for A=0.65 ,a'’= 0.1 5): V sin i = 3471±28 ms‘^ 
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Figure 5. Residuals of the Doppler data for HD 189733, computed in the vicinity of the transit epoch, corresponding to zero abscissa. 
Each raw of plots corresponds to a best fitting model labelled above the raw, together with the main best fitting parameters. See text 
for a detailed discussion. 


can prove a quick and practical workhorse, because most stars that are involved in planet search programmes are rather quiet, 
implying that their rotation should be relatively slow. _ _ _ 


Also, we do not consider the effect of macro-turbulence, which was considered e.g. by iHirano et al.l (1201 Ih : iBoue et al. 


(1201311 . and do not take into account differential rotation of the star. These and more subtle effects are left for future work. 

Regardless of all the remaining limitations, we believe that our model can be useful in practice, as it is fully analytic, 
requires nothing but the Doppler data, and can be applied without detailed knowledge of the spectrum reduction pipelines 
that depend on a particular practical case. This paper gives a comprehensive set of all necessary form ula, and we are going 
to release their implementation with the next version 3.0 of the PlanetPack package llBaluevll2013afl . 
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APPENDIX A: ROSSITER-MCLAUGHLIN ANOMALY FOR A SMALL PLANET, ARBITRARY 
ROTATION VELOCITY, AND MULTI-GAUSSIAN SPECTRA 

See Sect. [3] for the details of the approximation method. 


Al Cross-correlation with a predefined template 


Let us assume Gaussian approximation for all our spectra: 




(Al) 


Here the spectral lines positions u are the same for and but we admit that they may be slightly different from those 
used in the template, uj = u + Au. In this manner we model possible template imperfections. 

Using the expression and and approximating all slowly varying functions like Mq{s), and cy*(^) in the vicinity 

of each line by a constant, we can transform equation (Ireli to the following: 








(A2) 


If Am/ = 0 than this equality is satisfied automatically, and otherwise it sets a balancing requirement for Am/. For example, 
it is illegal if all Am/ only introduce a common systematic shift, because this would just result in a biasing effect on the RV 
absolute zero point, which does not affect relative RV measurements that we consider here. 

Various quantities appearing in (12.121) can be expressed analogously. Dropping the terms having relative magnitude 
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ff{Au^) and smaller, we obtain: 
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As we can see, the formulae for multiline spectra are significantly more complicated than for the single-line case considered 
in previous works. But before discussing them, let us consider the case of a single line. For A = 1 the formulae (IX3ll reduce to 
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(A4) 


This appears almost equivalent to the formula (15) by IIBoue et al.ll2013h . As they also fit the template via Ar, it becomes 
equal to At in their work. We obtain an additional fac tor of c/cj,, t he rat io of line intensities in the spectra of rotating star and 
stellar surface at rest. This ratio does not appear in Bone et al. 1 2013li . We believe this factor might be “lost” because they 
put an additional condition that and both should be pre-normalized, and consider them containing just a single line 
without even a continuum. This looks illegal, because the spectrum normalization mainly depends on its continuum, and not 
on the lines. Because line intensity a* may be different from a, the normalizations of and become different and cannot 
not be directly combined in Instead, it is better to consider unnormalized spectra treated as energy distributions, as we 
do in the present work. In this case we still do not need to take care of the continuum, but spectra normalizations become 
mutually consistent. 

The multiline approximation (TASll appears even more complicated. First, the summation over the lines in (IA3t should likely 
introduce additional broadening effect in comparison with the single-line formula (TAil) . Second, the multiline expression (IA31) 
contains terms depending on Am,-. This should introduce additional effect that depends on the quality of the template. This 
effect was not characterized previously, because it can be only revealed when working with the multiline model. Note that 
the functional shape of this template imperfection effect should be significantly different from the single-line formula (TAil) . 
Instead of the dependence on ip like G'{sp) ~ ipexp(—ip) (qualitatively), we should now deal with something like G”(ip) ~ 
(ip — l)exp(—ip). In fact, we cannot even guarantee that i = 0 for ip = 0 in this case: template imperfections introduce a bias. 

Unfortunately, this type of models is very difficult for practical use. It requires a comprehencive knowledge of deep 
internals of the spectra processing technique applied in the particular case. This is not available for authors who want to 
e.g. reanalyse some public Doppler data. Moreover, even when such a knowledge is available, the multiline model becomes 
mathematically complicated. Therefore, in this work we do not consider this type of approximations in more details. 

We do not give detailed expressions for the case in which the CCF is fitted by a Gaussian, described by the formula (12.1511 . 
Clearly, the hnal formulae for this case should be much more complicated than (I A3ll . Note that due to the template imperfec¬ 
tions, Au, appearing for multiline spectra, the term in (l2T5ll is non-zero even for symmetric lines and thus cannot 

be neglected. 


A2 Cross-correlation with an out-of-transit stellar spectrum or parametric modelling of the stellar 
spectrum (iodine cell technique) 

Now we should just substitute in place of in the formulae presented above. Formulae dMl) turn into 
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(A5) 


Here the template lines misplacements Am,- all vanish, because the new template coincides with Doppler anomaly can be 
expressed as follows: 
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Now the formula is more simple than (TASl) : the template imperfections are irrelevant, and the RM effect is not biased: Jp = 0 
implies f = 0. However, it still requires a very detailed knowledge of the stellar spectrum. 

For a single line, we obtain 
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'c,ypi + PiJ 2{p2 + p^)^ 

This basically agrees with Bone et al. 1 2013ll and Hirano et al. 1 2O10ll . But again, the summations in (IA 6 II should introduce 
an additional broadening effect in comparison with (IMt . 
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